<?php
/** 
* @package JAMA
*/

require_once "../Matrix.php";

/** 
* Example of use of Matrix Class, featuring magic squares. 
*/
class MagicSquareExample {

  /**
  * Generate magic square test matrix.
  * @param int n dimension of matrix
  */
  function magic($n) {
    
    // Odd order    
    
    if (($n % 2) == 1) {    
      $a = ($n+1)/2;
      $b = ($n+1);
      for ($j = 0; $j < $n; $j++) 
        for ($i = 0; $i < $n; $i++) 
          $M[$i][$j] = $n*(($i+$j+$a) % $n) + (($i+2*$j+$b) % $n) + 1;
    
    // Doubly Even Order
    
    } else if (($n % 4) == 0) {
      for ($j = 0; $j < $n; $j++) {
        for ($i = 0; $i < $n; $i++) {
          if ((($i+1)/2)%2 == (($j+1)/2)%2) 
            $M[$i][$j] = $n*$n-$n*$i-$j;
          else 
            $M[$i][$j] = $n*$i+$j+1;
        }
      }    
      
    // Singly Even Order      

    } else {                 

      $p = $n/2;  
      $k = ($n-2)/4;  
      $A = $this->magic($p);
      $M = array();
      for ($j = 0; $j < $p; $j++) {
        for ($i = 0; $i < $p; $i++) {
          $aij = $A->get($i,$j);
          $M[$i][$j]       = $aij;
          $M[$i][$j+$p]    = $aij + 2*$p*$p;
          $M[$i+$p][$j]    = $aij + 3*$p*$p;
          $M[$i+$p][$j+$p] = $aij + $p*$p;
        }
      }

      for ($i = 0; $i < $p; $i++) {
        for ($j = 0; $j < $k; $j++) { 
          $t = $M[$i][$j]; 
          $M[$i][$j] = $M[$i+$p][$j]; 
          $M[$i+$p][$j] = $t;
        }
        for ($j = $n-$k+1; $j < $n; $j++) {
          $t = $M[$i][$j]; 
          $M[$i][$j] = $M[$i+$p][$j]; 
          $M[$i+$p][$j] = $t;
        }
      }
      
      $t = $M[$k][0];  $M[$k][0]  = $M[$k+$p][0];  $M[$k+$p][0]  = $t;
      $t = $M[$k][$k]; $M[$k][$k] = $M[$k+$p][$k]; $M[$k+$p][$k] = $t;
    
    }
    
    return new Matrix($M);

  }

  /**
  * Simple function to replicate PHP 5 behaviour
  */
  function microtime_float() { 
    list($usec, $sec) = explode(" ", microtime()); 
    return ((float)$usec + (float)$sec); 
  } 

  /** 
  * Tests LU, QR, SVD and symmetric Eig decompositions.
  *
  *   n       = order of magic square.
  *   trace   = diagonal sum, should be the magic sum, (n^3 + n)/2.
  *   max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
  *   rank    = linear algebraic rank, should equal n if n is odd,
  *             be less than n if n is even.
  *   cond    = L_2 condition number, ratio of singular values.
  *   lu_res  = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
  *   qr_res  = test of QR factorization, norm1(Q*R-A)/(n*eps).
  */
  function main() {
    ?>
    <p>Test of Matrix Class, using magic squares.</p>
    <p>See MagicSquareExample.main() for an explanation.</p>          
    <table border='1' cellspacing='0' cellpadding='4'>
      <tr>
        <th>n</th>     
        <th>trace</th>       
        <th>max_eig</th>   
        <th>rank</th>        
        <th>cond</th>      
        <th>lu_res</th>           
        <th>qr_res</th>
      </tr>
      <?php
      $start_time = $this->microtime_float();       
      $eps = pow(2.0,-52.0);
      for ($n = 3; $n <= 6; $n++) {
        
        echo "<tr>";        
        
        echo "<td align='right'>$n</td>";        

        $M = $this->magic($n);
        $t = (int) $M->trace();                 
        
        echo "<td align='right'>$t</td>";  
        
        $O = $M->plus($M->transpose());
        $E = new EigenvalueDecomposition($O->times(0.5));       
        $d = $E->getRealEigenvalues();
        
        echo "<td align='right'>".$d[$n-1]."</td>";
        
        $r = $M->rank();
        
        echo "<td align='right'>".$r."</td>";        
        
        $c = $M->cond();
        
        if ($c < 1/$eps) 
          echo "<td align='right'>".sprintf("%.3f",$c)."</td>";                
        else           
          echo "<td align='right'>Inf</td>";                        

        $LU = new LUDecomposition($M);
        $L = $LU->getL();
        $U = $LU->getU();
        $p = $LU->getPivot();
        // Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1));                
        $S = $L->times($U);        
        $R = $S->minus($M->getMatrix($p,0,$n-1));                
        $res = $R->norm1()/($n*$eps);
        
        echo "<td align='right'>".sprintf("%.3f",$res)."</td>";                

        $QR = new QRDecomposition($M);
        $Q = $QR->getQ();
        $R = $QR->getR();
        $S = $Q->times($R);
        $R = $S->minus($M);
        $res = $R->norm1()/($n*$eps);
        
        echo "<td align='right'>".sprintf("%.3f",$res)."</td>";                         
        
        echo "</tr>";        
     
     }
     echo "<table>";     
     echo "<br />";

     $stop_time = $this->microtime_float();       
     $etime = $stop_time - $start_time;

     echo "<p>Elapsed time is ". sprintf("%.4f",$etime) ." seconds.</p>";

  }

}

$magic = new MagicSquareExample();
$magic->main();

?>
